

********************
*** Predict CPAP ***
********************

*** Sample Construction ***
{
*** MBSF - Part B 3 months, no HMO ***
use ".\temp\mbsf", clear
keep if year>=2009
keep bene_id year ZIP geo_id cba_round age gender ssdi enroll_Q* dual_Q*
reshape long enroll_Q dual_Q, i(bene_id age gender ssdi year geo_id cba_round) j(qtr)
keep if enroll_Q==3
gen yq = yq(year, qtr)
tempfile predict
save `predict'


*** Whether CPAP ***
use ".\statadata\dme_total_matched", clear
keep if hcpcs_cd=="E0601" & mod1=="RR" // CPAP
gen yq = qofd(clm_thru_dt)
collapse (sum) total_amt = line_alowd_chrg_amt q = line_srvc_cnt, by(bene_id yq)
gen CPAP = 1
merge 1:1 bene_id yq using `predict', nogen keep(2 3)
tempfile predict
save `predict'


*** OSA/SA Diagnosis ***
use ".\statadata\OSA_diag", clear
gen yq = qofd(date)
collapse (max) OSA = osa SA = sa, by(bene_id yq)
merge 1:1 bene_id yq using `predict', nogen keep(2 3)
foreach var in CPAP OSA SA {
    replace `var' = 0 if `var'==.
}

egen bene_idnum = group(bene_id)
xtset bene_idnum yq
gen enroll4 = 1 if (enroll_Q==3 & L.enroll_Q==3 & L2.enroll_Q==3 & L3.enroll_Q==3)
gen OSA4 = (OSA==1 | L.OSA==1 | L2.OSA==1 | L3.OSA==1) if enroll4==1
gen SA4 = (SA==1 | L.SA==1 | L2.SA==1 | L3.SA==1) if enroll4==1
tempfile predict
save `predict'
}


*** Logit ***
{
use `predict', clear
gen age_group = .
replace age_group = 0 if age<65
replace age_group = 1 if age>=65 & age<70
replace age_group = 2 if age>=70 & age<75
replace age_group = 3 if age>=75 & age<80
replace age_group = 4 if age>=80 & age<85
replace age_group = 5 if age>=85 & age<90
replace age_group = 6 if age>=90
gen dual = 1 if dual_Q==3 
replace dual = 0 if dual_Q==0
gen male = 1 if gender=="1"
replace male = 0 if gender=="2"

qui logit CPAP OSA4 SA4 male ib(1).age_group dual ssdi if year==2010 
	predict CPAP_hat_logit, pr
	outreg2 using $resdir\reg.dta, dta addstat(r2_p, e(r2_p)) replace

qui logit CPAP OSA4##dual SA4##dual male##dual ib(1).age_group##dual dual ssdi##dual if year==2010 
	predict CPAP_hat_logit2, pr
	outreg2 using $resdir\reg.dta, dta addstat(r2_p, e(r2_p))

preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA22_PredictCPAP") sheetreplace	
restore

tempfile predict
save `predict'
}



****************************
*** Event study by group ***
****************************
tempfile predict
save `predict'

foreach gg in 1 2 3 4 6  7 8 9 {
	use `predict', clear
	if `gg'==1 {
		gen group = dual
		local gg1 "Dual"
		local gg0 "Non-Dual"
	}
	if `gg'==2 {
		gen group = (male==0)
		local gg1 "Female"
		local gg0 "Male"
	}
	if `gg'==3 {
		gen group = (SA4==0) if SA!=.
		local gg1 "No SA Diagnosed"
		local gg0 "SA Diagnosed"
	}
	if `gg'==4 {
		gen group = (OSA4==0) if OSA!=.
		local gg1 "No OSA Diagnosed"
		local gg0 "OSA Diagnosed"
	}
	if `gg'==5 {
		gen temp = yq if CPAP==1
		bysort bene_id: egen fst_CPAP = min(temp)
		gen group = (fst_CPAP<yq)
		local gg1 "No CPAP Before"
		local gg0 "CPAP Before"
	}
	if `gg'==6 {
		bysort yq: egen temp = median(CPAP_hat_logit)
		gen group = (CPAP_hat_logit<temp) if CPAP_hat_logit!=.
		drop temp
		local gg1 "Below median CPAP propensity"
		local gg0 "Above median CPAP propensity"
	}
	if `gg'==7 {
	    keep if dual==1
		gen group = (SA4==0) if SA!=.
		local gg1 "No OSA Diagnosed"
		local gg0 "OSA Diagnosed"
	}
	if `gg'==8 {
	    keep if dual==1
		gen group = (OSA4==0) if OSA!=.
		local gg1 "No OSA Diagnosed"
		local gg0 "OSA Diagnosed"
	}
	if `gg'==9 {
	    keep if dual==1
		bysort yq: egen temp = median(CPAP_hat_logit)
		gen group = (CPAP_hat_logit<temp) if CPAP_hat_logit!=.
		drop temp
		local gg1 "Below median CPAP propensity"
		local gg0 "Above median CPAP propensity"
	}
	if `gg'==10 {
		bysort yq: egen temp = median(CPAP_hat_logit2)
		gen group = (CPAP_hat_logit2<temp) if CPAP_hat_logit2!=.
		drop temp
		local gg1 "Below median CPAP propensity"
		local gg0 "Above median CPAP propensity"
	}
	

	*** Regression sample construction ***
{
	gen pop = 1
	collapse (count) pop (sum) q total_amt, by(yq geo_id cba_round group)
	drop if yq==.
	drop if geo_id==""
	unique yq
	unique geo_id
	di 969*28*2
	tempfile temp
	save `temp'
	keep geo_id cba_round
	duplicates drop
	expand 28
	bysort geo_id: gen yq = yq(2008,4)+_n
	format yq %tq
	expand 2
	bysort geo_id yq: gen group = _n-1
	merge 1:1 geo_id yq group using `temp', nogen
	if `gg'==3 | `gg'==4 | `gg'==7 | `gg'==8 {
		keep if yq>=yq(2009,4)
	}
	if `gg'==6 | `gg'==9  | `gg'==10 {
		keep if yq>=yq(2011,1)
	}

	gen year = year(dofq(yq))
	merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen keepusing(pop_FFS_2010)

	replace q = 0 if q==.
	gen qX = q/pop*1000
	gen lg_qX = log(q)	
	gen price = total_amt/q
	gen lg_price = log(price)
	gen paymentX = total_amt/pop*1000
	gen lg_paymentX = log(total_amt)
}

	*** DID Regression by Subgroups ***
{
	gen dd = 0
	replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
	replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
	replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
	gen dd_group1 = dd*group
	gen dd_group0 = dd*(1-group)
	egen geo_idnum = group(geo_id)
	
	local rp
	if `gg'==1 | `gg'==7 {
		local rp replace
	}
	if `gg' < 7 {
	qui reghdfe lg_qX dd_group1 dd_group0 [aweight=pop_FFS_2010], nocons absorb(geo_idnum#group yq#group) vce(cluster geo_id)
		qui sum qX [aweight=pop_FFS_2010] if e(sample)==1 & group==1
		local mean1 = r(mean)
		qui sum qX [aweight=pop_FFS_2010] if e(sample)==1 & group==0
		local mean0 = r(mean)
		local diff = _b[dd_group1] - _b[dd_group0]
		qui test dd_group1 =  dd_group0
		local diff_p = r(p)
		outreg2 using $resdir\reg0.dta, dta addstat(subgroup,1, Mean of Dep Var (Group 1), `mean1', Mean of Dep Var (Group 2), `mean0', Coef Diff (Group 1 - 2), `diff', p-value, `diff_p') addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) `rp'
	local rp
	}
	if `gg' >=7 {
	qui reghdfe lg_qX dd_group1 dd_group0 [aweight=pop_FFS_2010], nocons absorb(geo_idnum#group yq#group) vce(cluster geo_id)
		qui sum qX [aweight=pop_FFS_2010] if e(sample)==1 & group==1
		local mean1 = r(mean)
		qui sum qX [aweight=pop_FFS_2010] if e(sample)==1 & group==0
		local mean0 = r(mean)
		local diff = _b[dd_group1] - _b[dd_group0]
		qui test dd_group1 =  dd_group0
		local diff_p = r(p)
		outreg2 using $resdir\reg2.dta, dta addstat(subgroup,1, Mean of Dep Var (Group 1), `mean1', Mean of Dep Var (Group 2), `mean0', Coef Diff (Group 1 - 2), `diff', p-value, `diff_p') addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) `rp'
	local rp
	}
}

	*** Event Study Plot ***
if `gg' < 7 {
	cap gen CBP = .
	replace CBP = yq(2011,1) if cba_round==1
	replace CBP = yq(2013,3) if cba_round==2 & CBP==.
	replace CBP = yq(2014,1) if cba_round==1 & CBP==.
	cap gen relative_qtr = yq-CBP
	tab relative_qtr, m
	drop if relative_qtr>9 & relative_qtr!=. 
	forvalues s = -20/9 {
		local q = `s'+21
		cap gen Q_`q' = (relative_qtr==`s')
		cap gen Q1_`q' = (relative_qtr==`s' & group==1)
		cap gen Q0_`q' = (relative_qtr==`s' & group==0)
	}

	local qtr_list Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 Q1_26 Q1_27 Q1_28 Q1_29 Q1_30  /// Q20 as baseline
		Q0_2 Q0_3 Q0_4 Q0_5 Q0_6 Q0_7 Q0_8 Q0_9 Q0_10 Q0_11 Q0_12 Q0_13 Q0_14 Q0_15 Q0_16 Q0_17 Q0_18 Q0_19 Q0_21 Q0_22 Q0_23 Q0_24 Q0_25 Q0_26 Q0_27 Q0_28 Q0_29 Q0_30 // Q20 (relative_qtr==-1 as baseline)

	foreach ver in lg_price lg_qX lg_paymentX {
		if "`ver'"=="lg_price" {
			local t "Log Price"
			local ylabel "ylabel(-0.8(0.2)0.2, valuelabel)"
		}
		if "`ver'"=="lg_qX" {
			local t "Log Quantity"
			local ylabel "ylabel(-0.3(0.1)0.2, valuelabel)"
		}
		if "`ver'"=="lg_paymentX" {
			local t "Log Spending"
			local ylabel "ylabel(-0.8(0.2)0.2, valuelabel)"
		}
		qui reghdfe `ver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(geo_idnum#group yq#group) vce(cluster geo_id)
		outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster replace
		
		preserve
		use $resdir\reg_dta.dta, clear
		gen i = 1 if v1!=""
		gen n = sum(i)
		bysort n: gen t = _n
		drop i
		reshape wide v1 v2, i(n) j(t)
		rename v11 var
		rename v21 coef
		rename v22 se
		rename v23 ci_low
		rename v24 ci_high
		gen relative_qtr = substr(var,4,.)
		gen group = substr(var,2,1)
		keep group relative_qtr coef se ci_low ci_high
		destring group relative_qtr coef se ci_low ci_high, replace force
		replace relative_qtr = relative_qtr-21
		drop if relative_qtr==.
		set obs 56
		replace relative_qtr = -1 if relative_qtr==.
		replace group = 1 in 55
		replace group = 0 in 56
		foreach s in se coef ci_low ci_high {
			replace `s' = 0 if relative_qtr==-1
		}
		tostring relative_qtr, gen(qtr)
		label define relative_qtr_v -20 "-20" -16 "-16" -12 "-12" -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8"
		label value relative_qtr relative_qtr_v
		sort group relative_qtr
		
		if (`gg'==1 | `gg'==3 | `gg'==4) & "`ver'"=="lg_qX" {
		twoway (connect coef relative_qtr if group==1, ms(o) mc(gs10) lc(gs10)) ///
			(rcap ci_low ci_high relative_qtr if group==1, lc(gs10)) ///
			(connect coef relative_qtr if group==0, ms(t) mc(black) lc(black) lp(dash)) ///
			(rcap ci_low ci_high relative_qtr if group==0, lc(black)), ///
			xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) `ylabel' ///
			legend(order(1 3) lab(1 "`gg1'") lab(3 "`gg0'") size(*0.8)) ///
			xtitle("Relative Qtr to CBP") title("`t'") graphregion(color(white))
		graph export $resdir\\Fig2_`ver'_`gg'.png, as(png) replace
		graph export $resdir\\Fig2_`ver'_`gg'.pdf, as(pdf) replace
		}
		if (`gg'==1 | `gg'==3 | `gg'==4) & ("`ver'"=="lg_price" | "`ver'"=="lg_paymentX") {
		twoway (connect coef relative_qtr if group==1, ms(o) mc(gs10) lc(gs10)) ///
			(rcap ci_low ci_high relative_qtr if group==1, lc(gs10)) ///
			(connect coef relative_qtr if group==0, ms(t) mc(black) lc(black) lp(dash)) ///
			(rcap ci_low ci_high relative_qtr if group==0, lc(black)), ///
			xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) `ylabel' ///
			legend(order(1 3) lab(1 "`gg1'") lab(3 "`gg0'") size(*0.8)) ///
			xtitle("Relative Qtr to CBP") title("`t'") graphregion(color(white))
		graph export $resdir\\FigA4_`ver'_`gg'.png, as(png) replace
		graph export $resdir\\FigA4_`ver'_`gg'.pdf, as(pdf) replace
		}
		
		if `gg'==6 | `gg'==10 {
			twoway (connect coef relative_qtr if group==1, ms(o) mc(gs10) lc(gs10)) ///
			(rcap ci_low ci_high relative_qtr if group==1, lc(gs10)) ///
			(connect coef relative_qtr if group==0, ms(t) mc(black) lc(black) lp(dash)) ///
			(rcap ci_low ci_high relative_qtr if group==0, lc(black)), ///
			xline(0, lwidth(vthin) lp(dash)) xlabel(-10(4)8, valuelabel) `ylabel' ///
			legend(order(1 3) lab(1 "`gg1'") lab(3 "`gg0'") size(*0.8)) ///
			xtitle("Relative Qtr to CBP") title("`t'") graphregion(color(white))
		graph export $resdir\\Fig3_`ver'_`gg'.png, as(png) replace
		graph export $resdir\\Fig3_`ver'_`gg'.pdf, as(pdf) replace
		}
		restore
	}

}

}

preserve
use $resdir\reg0_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table3_Subgroup") sheetreplace	
use $resdir\reg2_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA16_Subgroup") sheetreplace	
restore


* By Appropriateness 
{
use `predict', clear    

hist CPAP_hat_logit if SA4==1, percent bin(10) graphregion(color(white))
graph export $resdir\\FigRR_hist_CPAPhat_SA1.png, as(png) replace

hist CPAP_hat_logit if SA4==0, percent bin(10) graphregion(color(white))
graph export $resdir\\FigRR_hist_CPAPhat_SA0.png, as(png) replace
	
*xtile group = CPAP_hat_logit, nq(5)

gen group = .
sum CPAP_hat_logit if SA4==1, de
replace group = 4 if SA4==1 & CPAP_hat_logit>r(p50)
replace group = 3 if SA4==1 & CPAP_hat_logit<=r(p50)
sum CPAP_hat_logit if SA4==0, de
replace group = 2 if SA4==0 & CPAP_hat_logit>r(p50)
replace group = 1 if SA4==0 & CPAP_hat_logit<=r(p50)

gen pop = 1
collapse (count) pop (sum) q total_amt, by(yq geo_id cba_round group)
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28*4
tempfile temp
save `temp'
keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
expand 4
bysort geo_id yq: gen group = _n
merge 1:1 geo_id yq group using `temp', nogen
	
keep if yq>=yq(2011,1)

gen year = year(dofq(yq))
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen keepusing(pop_FFS_2010)

replace q = 0 if q==.
gen qX = q/pop*1000
gen lg_qX = log(q)	
gen price = total_amt/q
gen lg_price = log(price)

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
egen geo_idnum = group(geo_id)

local rp replace
forvalues i = 1/4 {
	qui reghdfe lg_qX dd [aweight=pop_FFS_2010] if group==`i', nocons absorb(geo_idnum yq) vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster `rp'
	local rp 
}

preserve 
use $resdir\reg_dta.dta, clear
gen i = 1 if v1!=""
gen n = sum(i)
bysort n: gen t = _n
drop i
keep if n==3
keep v2-v5
destring v2-v5, replace
xpose, clear varname
rename v1 coef
rename v2 se
rename v3 ci_low
rename v4 ci_high
gen group = substr(_varname,2,1)
destring group, replace
replace group = group - 1

twoway (scatter coef group, ms(o) mc(navy)) ///
	(rcap ci_low ci_high group, lc(navy)), ///
	yline(0, lwidth(vthin) lp(dash)) legend(off) ///
	xtitle("") xscale(r(0.7 4.3)) ///
	xlabel(1 `" "No diagnosis" "Low Propensity" "' 2 `" "No diagnosis" "High Propensity" "' ///
	3 `" "SA diagnosis" "Low Propensity" "' 4 `" "SA diagnosis" "High Propensity" "', labsize(small)) ///
	ytitle("Effect of CBP on CPAP Quantities") graphregion(color(white))
graph export $resdir\\FigRR_CPAP_lgqX_byApprop_SA_4.png, as(png) replace		
restore

}


* By Dual & Appropriateness 
{
use `predict', clear    

hist CPAP_hat_logit if dual==1, percent bin(10) graphregion(color(white))
graph export $resdir\\FigRR_hist_CPAPhat_dual.png, as(png) replace

hist CPAP_hat_logit if dual==0, percent bin(10) graphregion(color(white))
graph export $resdir\\FigRR_hist_CPAPhat_nondual.png, as(png) replace
	
*xtile group = CPAP_hat_logit, nq(5)

gen group = .
sum CPAP_hat_logit if dual==1, de
replace group = 4 if dual==1 & CPAP_hat_logit>r(p50)
replace group = 3 if dual==1 & CPAP_hat_logit<=r(p50)
sum CPAP_hat_logit if dual==0, de
replace group = 2 if dual==0 & CPAP_hat_logit>r(p50)
replace group = 1 if dual==0 & CPAP_hat_logit<=r(p50)

gen pop = 1
collapse (count) pop (sum) q total_amt, by(yq geo_id cba_round group)
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28*4
tempfile temp
save `temp'
keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
expand 4
bysort geo_id yq: gen group = _n
merge 1:1 geo_id yq group using `temp', nogen
	
keep if yq>=yq(2011,1)

gen year = year(dofq(yq))
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen keepusing(pop_FFS_2010)

replace q = 0 if q==.
gen qX = q/pop*1000
gen lg_qX = log(q)	
gen price = total_amt/q
gen lg_price = log(price)

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
egen geo_idnum = group(geo_id)

local rp replace
forvalues i = 1/4 {
	qui reghdfe lg_qX dd [aweight=pop_FFS_2010] if group==`i', nocons absorb(geo_idnum yq) vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster `rp'
	local rp 
}

preserve 
use $resdir\reg_dta.dta, clear
gen i = 1 if v1!=""
gen n = sum(i)
bysort n: gen t = _n
drop i
keep if n==3
keep v2-v5
destring v2-v5, replace
xpose, clear varname
rename v1 coef
rename v2 se
rename v3 ci_low
rename v4 ci_high
gen group = substr(_varname,2,1)
destring group, replace
replace group = group - 1

twoway (scatter coef group, ms(o) mc(navy)) ///
	(rcap ci_low ci_high group, lc(navy)), ///
	yline(0, lwidth(vthin) lp(dash)) legend(off) ///
	xtitle("") xscale(r(0.7 4.3)) ///
	xlabel(1 `" "Non Dual" "Low Propensity" "' 2 `" "Non Dual" "High Propensity" "' ///
	3 `" "Dual" "Low Propensity" "' 4 `" "Dual" "High Propensity" "', labsize(small)) ///
	ytitle("Effect of CBP on CPAP Quantities") graphregion(color(white))
graph export $resdir\\FigA7_CPAP_lgqX_byApprop_dual_4.png, as(png) replace		
restore

}



*************************
*** Marginal Consumer ***
*************************

* Regression
{
use `predict', clear
foreach ver of varlist CPAP_hat_logit* dual male age ssdi SA4 OSA4 {
	replace `ver' = . if CPAP==0
}
collapse (mean) CPAP CPAP_hat_logit* dual male age ssdi SA4 OSA4, by(yq geo_id cba_round)
drop if geo_id==""	
unique yq
unique geo_id
di 28*969

rename CPAP CPAP_rate
foreach ver of varlist CPAP_rate CPAP_hat_logit* dual male age ssdi SA4 OSA4 {
	cap gen lg_`ver' = log(`ver')
}
gen year = year(dofq(yq))
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen keepusing(pop_FFS_2010)

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
egen geo_idnum = group(geo_id)

local rp replace
foreach ver in dual male SA4 OSA4 {
	cap n qui ivreghdfe `ver' (lg_CPAP_rate = dd) [aweight=pop_FFS_2010], nocons absorb(geo_id yq) cluster(geo_id)
		qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in CPAP_hat_logit CPAP_hat_logit2 {
	cap n qui ivreghdfe `ver' (lg_CPAP_rate = dd) [aweight=pop_FFS_2010] if yq>=yq(2011,1), nocons absorb(geo_id yq) cluster(geo_id)
		qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Year FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean))
}	
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table5_MarginalConsumer") sheetreplace	
restore
}


* Event Study
{
gen CBP = .
replace CBP = yq(2011,1) if cba_round==1
replace CBP = yq(2013,3) if cba_round==2 & CBP==.
replace CBP = yq(2014,1) if cba_round==1 & CBP==.
		
gen relative_qtr = yq-CBP
tab relative_qtr, m
drop if (relative_qtr>8 & relative_qtr!=.) | relative_qtr<-8
forvalues i = -8/8 {
	local q = `i'+21
	gen Q`q' = (relative_qtr==`i')
}

replace CPAP_hat_logit = . if yq<yq(2011,1)

cap egen geo_idnum = group(geo_id)
local qtr_list Q13 Q14 Q15 Q16 Q17 Q18 Q19 /// Q20 as baseline
	 Q21 Q22 Q23 Q24 Q25 Q26 Q27 Q28 Q29  
	 
foreach ver in lg_CPAP_rate dual male SA4 OSA4 CPAP_hat_logit {
    if "`ver'"=="lg_CPAP_rate" {
		local tt "Log CPAP Rate"
	}
	if "`ver'"=="dual" {
		local tt "Dual"
	}
	if "`ver'"=="male" {
		local tt "Male"
	}
	if "`ver'"=="SA4" {
		local tt "Sleep Apnea (SA) Diagnosed"
	}
	if "`ver'"=="OSA4" {
		local tt "Obstructive Sleep Apnea (OSA) Diagnosed"
	}
	if "`ver'"=="CPAP_hat_logit" {
		local tt "CPAP Propensity"
	}
	qui reghdfe `ver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster replace
	
	preserve
	use $resdir\reg_dta.dta, clear
	gen i = 1 if v1!=""
	gen n = sum(i)
	bysort n: gen t = _n
	drop i
	reshape wide v1 v2, i(n) j(t)
	rename v11 var
	rename v21 coef
	rename v22 se
	rename v23 ci_low
	rename v24 ci_high
	keep var coef se ci_low ci_high
	gen relative_qtr = substr(var,2,.)
	keep relative_qtr coef se ci_low ci_high
	destring relative_qtr coef se ci_low ci_high, replace force
	replace relative_qtr = relative_qtr-21
	drop if relative_qtr==.
	set obs 28
	replace relative_qtr = -1 if relative_qtr==.
	foreach v in se coef ci_low ci_high {
		replace `v' = 0 if relative_qtr==-1
	}
	tostring relative_qtr, gen(qtr)
	label define relative_qtr_v -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8"
	label value relative_qtr relative_qtr_v
	sort relative_qtr
	scatter ci_high ci_low coef relative_qtr, connect(direct direct direct) ///
		lc(black black black) ms(none none none) lp(dot dot solid) ///
		xline(0, lwidth(vthin) lp(dash)) xlabel(-8(2)8, valuelabel) legend(off) ///
		xtitle("Relative Qtr to CBP") title("`tt'") graphregion(color(white))
	graph export $resdir\FigA8_`ver'.png, as(png) replace
	restore
}

}




*********************
*** Dual Analysis ***
*********************

* Sample Construction 
{
use ".\temp\sample_CPAP.dta", clear

bysort bene_id geo_id: egen temp = min(yq)
cap gen n_initiator = 1 if temp==yq
gen n_bene = 1 if q!=.

collapse (sum) n_bene n_initiator q total_amt n_claim, by(yq geo_id cba_round dual)
drop if dual==.
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28*2
tempfile temp
save `temp'
keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
expand 2
bysort geo_id yq: gen dual = _n-1
merge 1:1 geo_id yq dual using `temp', nogen

gen year = year(dofq(yq))
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen

replace q = 0 if q==.
replace n_bene = 0 if n_bene==.
replace n_initiator = 0 if n_initiator==.
replace total_amt = 0 if total_amt==.

gen qX = q/pop_FFSnondual*1000 if dual==0
replace qX = q/pop_FFSdual*1000 if dual==1
gen NbeneX = n_bene/pop_FFSnondual*1000 if dual==0
replace NbeneX = n_bene/pop_FFSdual*1000 if dual==1
gen NinitiatorX = n_initiator/pop_FFSnondual*1000 if dual==0
replace NinitiatorX = n_initiator/pop_FFSdual*1000 if dual==1
gen paymentX = total_amt/pop_FFSnondual*1000 if dual==0
replace paymentX = total_amt/pop_FFSdual*1000 if dual==1
	
gen lg_qX = log(qX)
gen lg_NbeneX = log(NbeneX)	
gen lg_NinitiatorX = log(NinitiatorX)	
gen lg_paymentX = log(paymentX)

gen price = total_amt/q
gen lg_price = log(price)

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
egen geo_idnum = group(geo_id)

cap gen nondual = (dual==0)
cap gen CBA_nondual = (cba_round>0 & nondual==1)
cap gen dd_nondual = dd*nondual
cap gen OOPC = 0.2*price if dual==0
replace OOPC = 0 if dual==1 & price!=.
}

* Regressions - Main
{
*** Table 6A: DDD ***
qui reghdfe OOPC dd_nondual dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id) 
qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) replace

foreach ver in qX NbeneX NinitiatorX {
	qui reghdfe lg_`ver' dd_nondual dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum `ver' if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) 
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table4_DualOLS") sheetreplace	
restore

*** Table 6B: IV ***
local rp replace
foreach ver in qX NbeneX NinitiatorX {
	qui ivreghdfe lg_`ver' (OOPC = dd_nondual) dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) cluster(geo_id)
	qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	local oopc_mean = r(mean)
	local elasticity = _b[OOPC]*`oopc_mean'
	local fstat = e(rkf)
	outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), `oopc_mean', Elasiticity, `elasticity', F stat, `fstat') `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table4_DualIV") sheetreplace	
restore

*** Table 6 RR: OLS ***
local rp replace
foreach ver in qX NbeneX NinitiatorX {
	qui reghdfe lg_`ver' OOPC dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum `ver' if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table4RR_OLS") sheetreplace	
restore
}

* State FEs
{
preserve
use ".\temp\ZIP_CBA_CBSA.dta", clear
keep geo_id CBA state state_name cba_round
drop if state=="" | geo_id==""
drop if geo_id=="CBA15" & state_name=="ME" // CBA = Boston, MA-NH only 1 zipcode
duplicates drop
gen R1 = 1 if cba_round==1
gen R2 = 1 if cba_round==2
bysort state: egen R1X = max(R1)
bysort state: egen R2X = max(R2)
drop R1 R2
expand 2
bysort geo_id state: gen dual = _n-1
expand 28
bysort geo_id state dual: gen yq = yq(2008,4)+_n
format yq %tq
levelsof state, local(levels)
foreach l of local levels {
	gen state_`l' = (state=="`l'")
	gen state_`l'_post_dual = 0
	replace state_`l'_post_dual = 1 if state=="`l'" & R1X==1 & yq>=yq(2011,1) & dual==1
	replace state_`l'_post_dual = 1 if state=="`l'" & R2X==1 & yq>=yq(2013,3) & dual==1
	
	gen state_`l'_post_nondual = 0
	replace state_`l'_post_nondual = 1 if state=="`l'" & R1X==1 & yq>=yq(2011,1) & dual==0
	replace state_`l'_post_nondual = 1 if state=="`l'" & R2X==1 & yq>=yq(2013,3) & dual==0
}
drop state_name
collapse (max) state_*, by(geo_id yq dual)
tempfile temp_state
save `temp_state'
restore
merge m:1 geo_id yq dual using `temp_state', nogen

* DDD
qui reg OOPC dd_nondual dd state_*_dual i.yq#i.dual i.geo_idnum#i.dual [aweight=pop_FFS_2010], vce(cluster geo_id)
qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
outreg2 using $resdir\reg.dta, dta keep(dd_nondual dd) addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, State-Dual-Post FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) replace

foreach ver in qX NbeneX NinitiatorX {
	qui reghdfe lg_`ver' dd_nondual dd state_*_dual [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum `ver' if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg.dta, dta keep(dd_nondual dd) addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, State-Dual-Post FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) 
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA19_DualOLS_state") sheetreplace	
restore

* IV
local rp replace
foreach ver in qX NbeneX NinitiatorX {
	qui ivreghdfe lg_`ver' (OOPC = dd_nondual) dd state_*_dual [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) cluster(geo_id)
	qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	local oopc_mean = r(mean)
	local elasticity = _b[OOPC]*`oopc_mean'
	local fstat = e(rkf)
	outreg2 using $resdir\reg.dta, keep(OOPC dd) dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, State-Dual-Post FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), `oopc_mean', Elasiticity, `elasticity', F stat, `fstat') `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA19_DualIV_state") sheetreplace	
restore

}

* State Medigap
{
preserve
use ".\statadata\mbsf", clear
keep if year==2015
gen ZIP = substr(zip_cd,1,5)
drop zip_cd
merge m:1 ZIP using ".\raw\ZIP_CBA_CBSA.dta", keep(1 3) nogen

drop if state=="" | geo_id==""
drop if bene_hmo_cvrage_tot_mons>0
keep geo_id CBA state state_name cba_round
gen n_pop = 1 
collapse (sum) n_pop, by(geo_id state_name)
tempfile temp_state
save `temp_state'

import excel "S:\group\mgduggan\CBP\raw\Medigap.xlsx", sheet("Sheet1") firstrow clear
rename State state_name
merge 1:m state_name using `temp_state', keep(3) nogen
collapse (mean) Share [aw=n_pop], by(geo_id)

sum Share, de
xtile share_2 = Share, nq(2)
xtile share_3 = Share, nq(3)
xtile share_4 = Share, nq(4)
replace share_4 = 2 if share_4==3 
replace share_4 = 3 if share_4==4
tempfile temp_state
save `temp_state'
restore
merge m:1 geo_id using `temp_state', nogen

local rp replace
forvalues i = 1/2 {
	qui reghdfe lg_qX dd_nondual dd [aweight=pop_FFS_2010] if share_2==`i', nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum qX if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg1.dta, dta keep(dd_nondual dd) addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean))  `rp'

	qui ivreghdfe lg_qX (OOPC = dd_nondual) dd [aweight=pop_FFS_2010] if share_2==`i', nocons absorb(geo_idnum#dual yq#dual) cluster(geo_id)
	qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	local oopc_mean = r(mean)
	local elasticity = _b[OOPC]*`oopc_mean'
	qui sum Share if e(sample) [aweight=pop_FFS_2010]
	local Medigap_share = r(mean)
	local fstat = e(rkf)
	outreg2 using $resdir\reg2.dta, keep(OOPC dd) dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), `oopc_mean', Elasiticity, `elasticity', Medigap Share, `Medigap_share', Medigap Median, `i', F stat, `fstat') `rp'
	local rp
}
forvalues i = 1/3 {
	qui reghdfe lg_qX dd_nondual dd [aweight=pop_FFS_2010] if share_3==`i', nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum qX if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg1.dta, dta keep(dd_nondual dd) addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) 

	qui ivreghdfe lg_qX (OOPC = dd_nondual) dd [aweight=pop_FFS_2010] if share_3==`i', nocons absorb(geo_idnum#dual yq#dual) cluster(geo_id)
	qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	local oopc_mean = r(mean)
	local elasticity = _b[OOPC]*`oopc_mean'
	qui sum Share if e(sample) [aweight=pop_FFS_2010]
	local Medigap_share = r(mean)
	local fstat = e(rkf)
	outreg2 using $resdir\reg2.dta, keep(OOPC dd) dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), `oopc_mean', Elasiticity, `elasticity', Medigap Share, `Medigap_share', Medigap Tercile, `i', F stat, `fstat')
}
preserve
use $resdir\reg1_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA17_DualOLS_state_medigap") sheetreplace	
use $resdir\reg2_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA17_DualIV_state_medigap") sheetreplace	
restore
}



*** Event Plot: OOPC ***
{
cap gen CBP = .
replace CBP = yq(2011,1) if cba_round==1
replace CBP = yq(2013,3) if cba_round==2 & CBP==.
replace CBP = yq(2014,1) if cba_round==1 & CBP==.
cap gen relative_qtr = yq-CBP
tab relative_qtr, m
drop if relative_qtr>9 & relative_qtr!=. 
forvalues s = -20/9 {
	local q = `s'+21
	cap gen Q_`q' = (relative_qtr==`s')
	cap gen Q1_`q' = (relative_qtr==`s' & dual==1)
	cap gen Q0_`q' = (relative_qtr==`s' & dual==0)
}

local qtr_list Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 Q1_26 Q1_27 Q1_28 Q1_29 Q1_30  /// Q20 as baseline
	Q0_2 Q0_3 Q0_4 Q0_5 Q0_6 Q0_7 Q0_8 Q0_9 Q0_10 Q0_11 Q0_12 Q0_13 Q0_14 Q0_15 Q0_16 Q0_17 Q0_18 Q0_19 Q0_21 Q0_22 Q0_23 Q0_24 Q0_25 Q0_26 Q0_27 Q0_28 Q0_29 Q0_30 // Q20 (relative_qtr==-1 as baseline)

foreach ver in OOPC {
	if "`ver'"=="OOPC" {
		local t "Out-of-Pocket Cost"
	}
	qui reghdfe `ver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) drop(CBA_dual) noparen noaster replace
	
	preserve
	use $resdir\reg_dta.dta, clear
	gen i = 1 if v1!=""
	gen n = sum(i)
	bysort n: gen t = _n
	drop i
	reshape wide v1 v2, i(n) j(t)
	rename v11 var
	rename v21 coef
	rename v22 se
	rename v23 ci_low
	rename v24 ci_high
	gen relative_qtr = substr(var,4,.)
	gen dual = substr(var,2,1)
	keep dual relative_qtr coef se ci_low ci_high
	destring dual relative_qtr coef se ci_low ci_high, replace force
	replace relative_qtr = relative_qtr-21
	drop if relative_qtr==.
	set obs 56
	replace relative_qtr = -1 if relative_qtr==.
	replace dual = 1 in 55
	replace dual = 0 in 56
	foreach s in se coef ci_low ci_high {
		replace `s' = 0 if relative_qtr==-1
	}
	tostring relative_qtr, gen(qtr)
	label define relative_qtr_v -20 "-20" -16 "-16" -12 "-12" -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8"
	label value relative_qtr relative_qtr_v
	sort dual relative_qtr
	twoway (connect coef relative_qtr if dual==1, ms(o) mc(navy) lc(navy)) ///
		(rcap ci_low ci_high relative_qtr if dual==1, lc(navy)) ///
		(connect coef relative_qtr if dual==0, ms(t) mc(maroon) lc(maroon) lp(dash)) ///
		(rcap ci_low ci_high relative_qtr if dual==0, lc(maroon)), ///
		xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) legend(order(1 3) lab(1 "Dual") lab(3 "Non-Dual")) ///
		xtitle("Relative Qtr to CBP") title("`t'") graphregion(color(white))
	graph export $resdir\\FigA5_`ver'.png, as(png) replace
	restore
}
}


**********************************************
*** Dual Analysis (state w/ full Medicaid) ***
**********************************************

* Sample Construction 
{
use ".\temp\sample_CPAP.dta", clear

bysort bene_id geo_id: egen temp = min(yq)
cap gen n_initiator = 1 if temp==yq
gen n_bene = 1 if q!=.

gen year = year(dofq(yq))
merge m:1 bene_id year using ".\temp\mbsf.dta", keepusing(ZIP) keep(1 3) nogen
merge m:1 ZIP using ".\temp\ZIP_CBA_CBSA.dta", keepusing(state_name) keep(1 3) nogen

gen state_sample = .
replace state_sample = 1 if inlist(state_name,"AR","HI","IA","ME","MO","MS") | inlist(state_name,"NE","OH","OK","SD","VT","WY") 
replace state_sample = 1 if inlist(state_name,"AL","AK","DE","IN","IA","MN","MT") | inlist(state_name,"NV","NM","NC","WV","WY")
keep if state_sample==1

collapse (sum) n_bene n_initiator q total_amt n_claim, by(yq geo_id cba_round dual)
drop if dual==.
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28*2
tempfile temp
save `temp'
keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
expand 2
bysort geo_id yq: gen dual = _n-1
merge 1:1 geo_id yq dual using `temp', nogen

gen year = year(dofq(yq))
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen

replace q = 0 if q==.
replace n_bene = 0 if n_bene==.
replace n_initiator = 0 if n_initiator==.
replace total_amt = 0 if total_amt==.

gen qX = q/pop_FFSnondual*1000 if dual==0
replace qX = q/pop_FFSdual*1000 if dual==1
gen NbeneX = n_bene/pop_FFSnondual*1000 if dual==0
replace NbeneX = n_bene/pop_FFSdual*1000 if dual==1
gen NinitiatorX = n_initiator/pop_FFSnondual*1000 if dual==0
replace NinitiatorX = n_initiator/pop_FFSdual*1000 if dual==1
gen paymentX = total_amt/pop_FFSnondual*1000 if dual==0
replace paymentX = total_amt/pop_FFSdual*1000 if dual==1
	
gen lg_qX = log(qX)
gen lg_NbeneX = log(NbeneX)	
gen lg_NinitiatorX = log(NinitiatorX)	
gen lg_paymentX = log(paymentX)

gen price = total_amt/q
gen lg_price = log(price)

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
egen geo_idnum = group(geo_id)

cap gen nondual = (dual==0)
cap gen CBA_nondual = (cba_round>0 & nondual==1)
cap gen dd_nondual = dd*nondual
cap gen OOPC = 0.2*price if dual==0
replace OOPC = 0 if dual==1 & price!=.
}

* Regressions - Main
{
*** Table 6A: DDD ***
qui reghdfe OOPC dd_nondual dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id) 
qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) replace

foreach ver in qX NbeneX NinitiatorX {
	qui reghdfe lg_`ver' dd_nondual dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	qui sum `ver' if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), r(mean)) 
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA20_DualOLS_fullMCD") sheetreplace	
restore

*** Table 6B: IV ***
local rp replace
foreach ver in qX NbeneX NinitiatorX {
	qui ivreghdfe lg_`ver' (OOPC = dd_nondual) dd [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) cluster(geo_id)
	qui sum OOPC if e(sample) & year==2010 & dual==0 [aweight=pop_FFS_2010]
	local oopc_mean = r(mean)
	local elasticity = _b[OOPC]*`oopc_mean'
	local fstat = e(rkf)
	outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr-Dual FE, X, GEO-Dual FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var (Non-dual Pre), `oopc_mean', Elasiticity, `elasticity', F stat, `fstat') `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA20_DualIV_fullMCD") sheetreplace	
restore
}

